cubble::climate_flat

## # A tibble: 1,830 × 10
##   id         lat  long  elev name  wmo_id date        prcp  tmax  tmin
##   <chr>    <dbl> <dbl> <dbl> <chr>  <dbl> <date>     <dbl> <dbl> <dbl>
## 1 ASN0000… -31.9  116.  15.4 pert…  94610 2020-01-01     0  31.9  15.3
## 2 ASN0000… -31.9  116.  15.4 pert…  94610 2020-01-02     0  24.9  16.4
## 3 ASN0000… -31.9  116.  15.4 pert…  94610 2020-01-03     6  23.2  13  
## 4 ASN0000… -31.9  116.  15.4 pert…  94610 2020-01-04     0  28.4  12.4
## 5 ASN0000… -31.9  116.  15.4 pert…  94610 2020-01-05     0  35.3  11.6
## # … with 1,825 more rows

cubble_nested <- climate_flat |>
  as_cubble(key = id, index = date, coords = c(long, lat))
cubble_nested

## # cubble:   id [5]: nested form
## # bbox:     [115.97, -32.94, 133.55, -12.42]
## # temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
##   id            lat  long  elev name           wmo_id ts      
##   <chr>       <dbl> <dbl> <dbl> <chr>           <dbl> <list>  
## 1 ASN00009021 -31.9  116.  15.4 perth airport   94610 <tibble>
## 2 ASN00010311 -31.9  117. 179   york            94623 <tibble>
## 3 ASN00010614 -32.9  117. 338   narrogin        94627 <tibble>
## 4 ASN00014015 -12.4  131.  30.4 darwin airport  94120 <tibble>
## 5 ASN00015131 -17.6  134. 220   elliott         94236 <tibble>

cubble_nested |>
  mutate(rain_day = sum(ts$prcp != 0))

## # cubble:   id [5]: nested form
## # bbox:     [115.97, -32.94, 133.55, -12.42]
## # temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
##   id            lat  long  elev name          wmo_id ts       rain_day
##   <chr>       <dbl> <dbl> <dbl> <chr>          <dbl> <list>      <int>
## 1 ASN00009021 -31.9  116.  15.4 perth airport  94610 <tibble>      104
## 2 ASN00010311 -31.9  117. 179   york           94623 <tibble>       89
## 3 ASN00010614 -32.9  117. 338   narrogin       94627 <tibble>       90
## 4 ASN00014015 -12.4  131.  30.4 darwin airpo…  94120 <tibble>      106
## 5 ASN00015131 -17.6  134. 220   elliott        94236 <tibble>       63

cubble_long <- cubble_nested |> face_temporal()
cubble_long

## # cubble:  date, id [5]: long form
## # bbox:    [115.97, -32.94, 133.55, -12.42]
## # spatial: lat [dbl], long [dbl], elev [dbl], name [chr], wmo_id
## #   [dbl]
##   id          date        prcp  tmax  tmin
##   <chr>       <date>     <dbl> <dbl> <dbl>
## 1 ASN00009021 2020-01-01     0  31.9  15.3
## 2 ASN00009021 2020-01-02     0  24.9  16.4
## 3 ASN00009021 2020-01-03     6  23.2  13  
## 4 ASN00009021 2020-01-04     0  28.4  12.4
## 5 ASN00009021 2020-01-05     0  35.3  11.6
## # … with 1,825 more rows

cubble_back <- cubble_long |> face_spatial()
cubble_back

## # cubble:   id [5]: nested form
## # bbox:     [115.97, -32.94, 133.55, -12.42]
## # temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
##   id            lat  long  elev name           wmo_id ts      
##   <chr>       <dbl> <dbl> <dbl> <chr>           <dbl> <list>  
## 1 ASN00009021 -31.9  116.  15.4 perth airport   94610 <tibble>
## 2 ASN00010311 -31.9  117. 179   york            94623 <tibble>
## 3 ASN00010614 -32.9  117. 338   narrogin        94627 <tibble>
## 4 ASN00014015 -12.4  131.  30.4 darwin airport  94120 <tibble>
## 5 ASN00015131 -17.6  134. 220   elliott         94236 <tibble>

identical(cubble_nested, cubble_back)

## [1] TRUE

cubble_unfold <- cubble_long |> unfold(long, lat)
cubble_unfold

## # cubble:  date, id [5]: long form
## # bbox:    [115.97, -32.94, 133.55, -12.42]
## # spatial: lat [dbl], long [dbl], elev [dbl], name [chr], wmo_id
## #   [dbl]
##   id          date        prcp  tmax  tmin  long   lat
##   <chr>       <date>     <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ASN00009021 2020-01-01     0  31.9  15.3  116. -31.9
## 2 ASN00009021 2020-01-02     0  24.9  16.4  116. -31.9
## 3 ASN00009021 2020-01-03     6  23.2  13    116. -31.9
## 4 ASN00009021 2020-01-04     0  28.4  12.4  116. -31.9
## 5 ASN00009021 2020-01-05     0  35.3  11.6  116. -31.9
## # … with 1,825 more rows

climate_flat_ts <- climate_flat |> 
  tsibble::as_tsibble(key = id, index = date)
climate_flat_cb <-  climate_flat_ts |> 
  cubble::as_cubble(coords = c(long, lat))
climate_flat_cb

## # cubble:   id [5]: nested form
## # bbox:     [115.97, -32.94, 133.55, -12.42]
## # temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
##   id            lat  long  elev name           wmo_id ts      
##   <chr>       <dbl> <dbl> <dbl> <chr>           <dbl> <list>  
## 1 ASN00009021 -31.9  116.  15.4 perth airport   94610 <tbl_ts>
## 2 ASN00010311 -31.9  117. 179   york            94623 <tbl_ts>
## 3 ASN00010614 -32.9  117. 338   narrogin        94627 <tbl_ts>
## 4 ASN00014015 -12.4  131.  30.4 darwin airport  94120 <tbl_ts>
## 5 ASN00015131 -17.6  134. 220   elliott         94236 <tbl_ts>

climate_flat_cb |> 
  mutate(
    fabletools::features(
      ts, tmax, list(tmax_mean = mean, tmax_var = var)
      )
  )

## # cubble:   id [5]: nested form
## # bbox:     [115.97, -32.94, 133.55, -12.42]
## # temporal: date [date], prcp [dbl], tmax [dbl], tmin [dbl]
##   id          lat  long  elev name  wmo_id ts       tmax_mean tmax_var
##   <chr>     <dbl> <dbl> <dbl> <chr>  <dbl> <list>       <dbl>    <dbl>
## 1 ASN00009… -31.9  116.  15.4 pert…  94610 <tbl_ts>      25.7    38.6 
## 2 ASN00010… -31.9  117. 179   york   94623 <tbl_ts>      26.2    51.1 
## 3 ASN00010… -32.9  117. 338   narr…  94627 <tbl_ts>      23.7    45.4 
## 4 ASN00014… -12.4  131.  30.4 darw…  94120 <tbl_ts>      33.1     3.02
## 5 ASN00015… -17.6  134. 220   elli…  94236 <tbl_ts>      34.6    24.7

## library(tidyverse)
## #https://www.coronavirus.vic.gov.au/victorian-coronavirus-covid-19-data
## raw <- read_csv(here::here("data/NCOV_COVID_Cases_by_LGA_Source_20220324.csv"),
##                   col_names = c("date", "postcode", "source", "lga"), skip = 1)
## 
## covid <- raw |>
##   count(date, lga, source, postcode) |>
##   filter(date >= as.Date("2022-01-01")) |>
##   filter(source == "Contact with a confirmed case") |>
##   group_by(date, lga, source) |>
##   count() |>
##   ungroup(date) |>
##   mutate(roll_mean = zoo::rollmean(n, k = 14, fill = NA)) |>
##   tsibble::as_tsibble(key = lga, index = date) |>
##   tsibble::fill_gaps(.full = TRUE)
## save(covid, file = here::here("data/covid.rda"))

load(here::here("data/covid.rda"))

lga <- strayr::read_absmap("lga2018") |>
  rename(lga = lga_name_2018) |>
  filter(state_name_2016 == "Victoria") 

covid |> head(5)

## # A tsibble: 5 x 5 [1D]
## # Key:       lga [1]
## # Groups:    lga, source [1]
##   date       lga        source                            n roll_mean
##   <date>     <chr>      <chr>                         <int>     <dbl>
## 1 2022-01-01 Alpine (S) Contact with a confirmed case     1        NA
## 2 2022-01-02 Alpine (S) Contact with a confirmed case     2        NA
## 3 2022-01-03 Alpine (S) Contact with a confirmed case     4        NA
## 4 2022-01-04 Alpine (S) Contact with a confirmed case     4        NA
## 5 2022-01-05 Alpine (S) Contact with a confirmed case     2        NA

lga |> head(5)

## Simple feature collection with 5 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 142.3535 ymin: -38.67876 xmax: 147.3909 ymax: -36.39269
## Geodetic CRS:  WGS 84
##     lga_code_2018            lga state_code_2016 state_name_2016
## 132         20110     Alpine (S)               2        Victoria
## 133         20260    Ararat (RC)               2        Victoria
## 134         20570   Ballarat (C)               2        Victoria
## 135         20660    Banyule (C)               2        Victoria
## 136         20740 Bass Coast (S)               2        Victoria
##     areasqkm_2018 cent_long  cent_lat                       geometry
## 132     4788.1568  146.9742 -36.85357 MULTIPOLYGON (((146.7258 -3...
## 133     4211.1171  142.8432 -37.47271 MULTIPOLYGON (((143.1807 -3...
## 134      739.0321  143.7815 -37.49286 MULTIPOLYGON (((143.6622 -3...
## 135       62.5402  145.0851 -37.73043 MULTIPOLYGON (((145.1357 -3...
## 136      865.8095  145.5581 -38.50730 MULTIPOLYGON (((145.5207 -3...

cb <- as_cubble(
  list(spatial = lga, temporal = covid),
  key = lga, index = date, coords = c(cent_long, cent_lat)
  )

## ! Some sites in the temporal table don't have corresponding spatial information

## ! Some sites in the spatial table don't have corresponding temporal information

## ! Use argument `output = "unmatch"` to check on the unmatched key

pair <- as_cubble(
  list(spatial = lga, temporal = covid),
  key = lga, index = date, coords = c(cent_long, cent_lat),
  output = "unmatch"
  )

pair

## $paired
## # A tibble: 2 × 2
##   spatial             temporal    
##   <chr>               <chr>       
## 1 Kingston (C) (Vic.) Kingston (C)
## 2 Latrobe (C) (Vic.)  Latrobe (C) 
## 
## $others
## $others$temporal
## [1] "Interstate" "Overseas"   "Unknown"   
## 
## $others$spatial
## [1] "No usual address (Vic.)"               
## [2] "Migratory - Offshore - Shipping (Vic.)"

lga <- lga |>
  mutate(lga = ifelse(lga == "Kingston (C) (Vic.)", "Kingston (C)", lga),
         lga = ifelse(lga == "Latrobe (C) (Vic.)", "Latrobe (C)", lga)) |>
  filter(!lga %in% pair$others$spatial)

covid <- covid |> filter(!lga %in% pair$others$temporal)

cb <- as_cubble(data = list(spatial = lga, temporal = covid),
                key = lga, index = date, coords = c(cent_long, cent_lat))

load(here::here("data/historical_tmax.rda"))

tmax <- historical_tmax |>
  filter(between(stringr::str_sub(id, 7, 8), 46, 90))

tmax <- tmax |>
  face_temporal() |> 
  group_by(month = lubridate::month(date),
         group = as.factor(
           ifelse(lubridate::year(date) > 2015,
           "2016 ~ 2020", "1971 ~ 1975"))) |>
  summarise(tmax = mean(tmax, na.rm = TRUE))

tmax |> 
  face_spatial() |> 
  mutate(n = nrow(ts)) |>
  arrange(n) |>
  pull(n) |> 
  head(10)

##  [1] 12 12 12 13 19 19 20 24 24 24

tmax <- tmax |> 
  face_spatial() |> 
  filter(nrow(ts) == 24) |>
  face_temporal()

cobar <- tmax |> filter(id == "ASN00048027") |>
  ggplot(aes(x = month, y = tmax, color = group)) +
  geom_line() +
  scale_color_brewer(palette = "Dark2", guide = "none") +
  scale_x_continuous(
    breaks = seq(1, 12, 1), 
    labels = c("J", "F", "M", "A",
               "M", "J", "J", "A",
               "S", "O", "N", "D")
    ) +
  labs(x = "", y  = "Temp (C)") +
  theme_bw() + 
  theme(aspect.ratio = 0.3)

tmax <- tmax |> unfold(long, lat)
box_df <- tmax |> face_spatial() |> filter(id == "ASN00048027")
nsw_vic <- ozmaps::abs_ste |> filter(NAME %in% c("Victoria","New South Wales"))

p1 <- ggplot() + 
  geom_sf(data = nsw_vic, fill = "grey95", 
          color = "white", size=2) + 
  geom_glyph(
    data = tmax, 
    aes(x_major = long, x_minor = month, 
        y_major = lat, y_minor = tmax,
        group = interaction(id, group), color = group),
    width = 1.2, height = 0.8) +
    geom_rect(
      data = box_df,
      aes(xmin = long - 0.6, xmax = long + 0.6,
          ymin = lat - 0.12, ymax = lat + 0.35),
      fill = "transparent", color = "black") +
  scale_color_brewer("", palette = "Dark2") + 
  coord_sf(xlim = c(141, 154), ylim = c(-39, -28.5)) + 
  ggthemes::theme_map() +
  theme(legend.position = "bottom", legend.text = element_text(size = 10)) + 
  guides(color = guide_legend(override.aes = list(size=2)))

(p1 | cobar)  + 
  patchwork::plot_layout(width = c(1,1)) + 
  plot_annotation(tag_levels = 'a')

A glyph map of the monthly maximum average temperature for weather stations in Victoria and New South Wales (NSW) for the periods (1971-1975, 2016-2020). The corresponding average time series for the cobar station are display on the top left corner. From the glyph map we can observe that the monthly trend is similar for all locations (low in the winter, high in the summer), and small increased temperatures, particularly in late summer can be seen at most stations in NSW.

load(here::here("data/climate_full.rda"))

climate <- climate_full |>
  filter(between(stringr::str_sub(id, 7, 8), 76, 90)) |>
  face_temporal() |>
  filter(lubridate::year(date) == 2020) |>
  face_spatial() |>
  mutate(type = "climate")

river <- river |> mutate(type = "river")

res <- match_sites(
  river, climate,
  temporal_by = c("Water_course_level" = "prcp"),
  temporal_independent = "prcp",  
  temporal_n_highest = 30,
  temporal_min_match = 15, 
)

res |> print()

## # cubble:   id [8]: nested form
## # bbox:     [144.52, -37.73, 146.06, -36.55]
## # temporal: date [date], matched_var [dbl]
##   id          name        lat  long type   dist group ts       n_match
##   <chr>       <chr>     <dbl> <dbl> <chr> <dbl> <int> <list>     <int>
## 1 405234      SEVEN CR… -36.9  146. river  6.15     5 <tibble>      21
## 2 404207      HOLLAND … -36.6  146. river  8.54    10 <tibble>      21
## 3 ASN00082042 strathbo… -36.8  146. clim…  6.15     5 <tibble>      21
## 4 ASN00082170 benalla … -36.6  146. clim…  8.54    10 <tibble>      21
## 5 230200      MARIBYRN… -37.7  145. river  6.17     6 <tibble>      19
## # … with 3 more rows

vic_map <- ozmaps::abs_ste |> 
  filter(NAME == "Victoria") |> 
  rmapshaper::ms_simplify()  

p1 <-ggplot() + 
  geom_sf(data = vic_map, fill = "grey95", color = "white") + 
  geom_point(data = dplyr::bind_rows(river, climate), 
             aes(x = long, y = lat, color = type), alpha = 0.2, fill = 0.2) +
  geom_point(data = res, aes(x = long, y = lat, color = type)) +
  ggrepel::geom_label_repel(
    data = res |> filter(type == "river"), aes(x = long, y = lat, label = group)) +
  scale_color_brewer(palette = "Dark2")  + 
  ggthemes::theme_map() +
  ggplot2::theme(legend.position = "bottom",
                 legend.text = element_text(size = 15),
                 legend.title = element_text(size = 15)) +
  ggplot2::labs(x = "Longitude", y = "Latitude") + 
  guides(color = guide_legend(override.aes = list(size=5)))

res_long <- res |>
  face_temporal(ts) |>
  unfold(group, type) |>
  rename(prcp = matched_var) |> 
  mutate(prcp = (prcp - min(prcp, na.rm = TRUE))/ (max(prcp, na.rm = TRUE) - min(prcp, na.rm = TRUE))) 

p2 <- res_long |> 
  ggplot(aes(x = date, y = prcp, color = type, group = id)) +
  geom_line(alpha = 0.8) +
  facet_wrap(vars(group)) +
  scale_color_brewer(palette = "Dark2", guide = "none") +
  theme_bw() +
  scale_x_date(date_labels = "%b") + 
  labs(x = "Date", y = "Precipitation/ water level")

(p1 | p2) + 
  patchwork::plot_layout(guides = "collect") + 
  plot_annotation(tag_levels = 'a')&
  theme(legend.position = "bottom") 

Weather stations and river gauges with matched pairs labelled on the map (a) and plotted across time (b). Precipitation and water level have been standardised between 0 and 1 to be displayed on the same scale. The water level reflects the increase in precipitation. The numbers (1, 5, 6, 10) indicate the group index derived from spatial matching, only those that were selectd by temporal matching are shown here.

raw <- ncdf4::nc_open(here::here("data/era5-pressure.nc"))
dt <- as_cubble(
  raw, vars = c("q", "z"),
  long_range = seq(-180, 180, 1), lat_range = seq(-88, 88, 1))

## # prepare data/era5-pressure.nc
## library(ecmwfr)
## my_request <- list(
##   product_type = "reanalysis",
##   format = "netcdf",
##   pressure_level = "10",
##   year = "2002",
##   month = c("09", "10"),
##   day = c("04", "22", "26", "30"),
##   time = "00:00",
##   area = c(-5, -180, -90, 180),
##   variable = c("specific_humidity", 'geopotential'),
##   dataset_short_name = "reanalysis-era5-pressure-levels",
##   target = "era5-pressure.nc"
## )
## #
## # # here you need to create an account on Climate Data Store
## # # at https://cds.climate.copernicus.eu/cdsapp#!/home) to obtain the user ID and API key
## # # then set up the with ecmwfr::wf_set_key(user, key, service = "cds")
## # # the user below is the same `user` in the setup line
## wf_request(
##   user = "113339",
##   request = my_request,
##   transfer = TRUE,
##   path = here::here("data/"))

date <- c("2002-09-22", "2002-09-26", "2002-09-30", "2002-10-04") |> as.Date()
res <- dt |> 
  face_temporal() |> 
  filter(lubridate::date(time) %in% date) |>
  unfold(long, lat) |> 
  mutate(q =  q* 10^6)

box = c(xmin = -180, ymin = -90, xmax = 180, ymax = -5)
sf::sf_use_s2(FALSE)
country <- rnaturalearth::ne_coastline("small", returnclass = "sf") |> 
  sf::st_crop(box) |> 
  sf::st_cast("MULTILINESTRING")

res |> 
  ggplot() +
  # q for specific humidity
  geom_point(aes(x = long, y = lat, color = q), size = 2) + 
  #geom_tile(aes(x = long, y = lat, fill = q), color = "transparent") +
  # z for geopotential
  geom_contour(data = res, aes(x = long, y = lat, z = z),
               color = "grey20", binwidth = 4000, size = 0.5) +
  geom_sf(data = country , alpha = 0.5, fill = "transparent", color = "lightgreen") +
  coord_sf(
    crs = "+proj=stere +lat_0=-90 +lon_0=-180 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ", default_crs = sf::st_crs(4326), label_graticule = "EW") + 
  theme_bw() + 
  theme(legend.position = "bottom",
        axis.ticks = element_blank(),
        panel.border = element_blank(),
        axis.text = element_blank(), 
        panel.grid = element_blank()) + 
  facet_wrap(vars(as.Date(time)), nrow = 1, shrink = TRUE)+  
  colorspace::scale_color_continuous_sequential("Purple-Yellow", name = "Specific humidity") + 
  theme(legend.text = element_text(size = 13),
        legend.title = element_text(size = 13)) + 
  labs(x ="", y = "")

A reproduction of the second row (ERA5 data) of Figure 19 in Hersbach et al (2020) to illustrate the break-up of sourthern polar vortex in late September and early October 2002. The polar vortex, signalled by the high speicfic humidity, splits into two on 2002-09-26 and furthers split into four on 2002-10-04.

## # This chunk generates the interactive graphic and save it as figures/linking.html, where the three screencasts in the figure comes from.
## library("crosstalk")
## library("plotly")
## 
## clean <- climate_full |>
##   filter(!id %in% c("ASN00067033", "ASN00072091", "ASN00059040",
##                     "ASN00097053", "ASN00040856", "ASN00015548")) |>
##   face_temporal() |>
##   mutate(month = lubridate::month(date)) |>
##   group_by(month) |>
##   summarise(
##     tmax = mean(tmax, na.rm = TRUE),
##     tmin = mean(tmin, na.rm = TRUE),
##     diff = mean(tmax - tmin, na.rm = TRUE)
##     ) |>
##   ungroup(month) |>
##   mutate(dummy_date = as.Date(glue::glue("2021-{month}-01"))) |>
##   face_spatial() |>
##   mutate(temp_diff_var = ts |> pull(diff) |> var(na.rm = TRUE))
## 
## nested <- clean |>
##   SharedData$new(~id, group = "cubble")
## 
## long <- clean |>
##   face_temporal() |>
##   unfold(temp_diff_var) |>
##   arrange(temp_diff_var) |>
##   SharedData$new(~id, group = "cubble")
## 
## state_map <- rmapshaper::ms_simplify(ozmaps::abs_ste, keep = 2e-3)
## 
## p1 <- nested |>
##   ggplot() +
##   geom_sf(data = state_map, fill = "grey95", color = "white") +
##   geom_point(aes(x = long, y = lat,color = temp_diff_var, label = name)) +
##   colorspace::scale_color_continuous_sequential(
##     "Rocket",  n_interp = 7, cmax = 90, rev = TRUE,
##     c2 = 40, l2= 85, c1 = 20, l1 = 30, name = "Var. temp. diff.") +
##   labs(x = "Longitude", y = "Latitude") +
##   coord_sf() +
##   theme_bw() +
##   theme(
##     panel.grid.major = element_blank(),
##     legend.position = "none"
##     )
## 
## p2 <- long |>
##   ggplot(aes(x = dummy_date, group = id,
##              fill = temp_diff_var, color = temp_diff_var)) +
##   geom_ribbon(aes(ymin = tmin, ymax = tmax), size = 0.1, alpha = 0.3) +
##   geom_line(aes(y = tmin), alpha = 0.1) +
##   geom_line(aes(y = tmax), alpha = 0.1) +
##   geom_point(aes(y = tmax), size = 0.1) +
##   geom_point(aes(y = tmin), size = 0.1) +
##   colorspace::scale_fill_continuous_sequential(
##     "Rocket",  n_interp = 7, cmax = 90, rev = TRUE,
##     c2 = 40, l2= 85, c1 = 20, l1 = 30, name = "Var. temp. diff.") +
##   colorspace::scale_colour_continuous_sequential(
##     "Rocket",  n_interp = 7, cmax = 90, rev = TRUE,
##     c2 = 40, l2= 85, c1 = 20, l1 = 30, name = "Var. temp. diff.") +
##   scale_x_date(date_labels = "%b", date_breaks = "1 month") +
##   labs(x = "Month", y = "Temperature") +
##   theme_bw() +
##   theme(
##     panel.grid.major = element_blank(),
##     legend.position = "bottom"
##     )
## 
## out <- bscols(
##   ggplotly(p1, width = 900, height = 700, tooltip = "label") |>
##     highlight(
##       on = "plotly_selected",
##       off = "plotly_deselect",
##       opacityDim = 0.1),
##   ggplotly(p2, width = 1200, height = 700, tooltip = "label") |>
##     highlight(
##       on = "plotly_selected",
##       off = "plotly_deselect",
##       opacityDim = 0.012),
##   widths = c(5, 5)
## )
## 
## htmltools::save_html(out, file = "figures/linking.html")
## 

library("leaflet")
library("leafpop")

## # The kmeans clustering in station_nested data can takes some time to run,
## # This chunk runs the algorithm and save the result as staiton_nested for further use
## coords <- cbind(prcp_aus$long, prcp_aus$lat)
## dist_raw <- geosphere::distm(coords, coords)
## 
## station_long <- cubble::prcp_aus |>
##   face_temporal(ts) |>
##   mutate(wk = lubridate::week(date)) |>
##   group_by(wk) |>
##   summarise(prcp = sum(prcp, na.rm = TRUE)) |>
##   unfold(cluster)
## 
## set.seed(123)
## station_nested <- station_long |>
##   face_spatial() |>
##   strip_rowwise() |>
##   mutate(cluster = kmeans(dist_raw,centers = 20, nstart = 500)$cluster)
## save(station_nested, file = here::here("data/station_nested.rda"))

load(here::here("data/station_nested.rda"))
cluster_nested <- station_nested |>
  switch_key(cluster) |> 
  get_centroid()

cluster_long <- cluster_nested |> 
  face_temporal() |> 
  group_by(wk) |>
  summarise(prcp = mean(prcp, na.rm = TRUE)) |> 
  unfold(cent_long, cent_lat)

state_map <- rmapshaper::ms_simplify(ozmaps::abs_ste, keep = 2e-3)

ggplot_smooth <- cluster_long |> 
  ggplot() +
  geom_smooth(aes(x = wk, y = prcp, group = cluster), span = 0.4) 

smoother <- layer_data(ggplot_smooth) |> 
  left_join(cluster_long |> select(cluster, cent_long, cent_lat), by = c("group" = "cluster"))

p1 <- ggplot(data = smoother, 
             aes(x_minor = x, y_minor = y, x_major = cent_long, y_major = cent_lat)) + 
  geom_sf(data = state_map, inherit.aes = FALSE,
          color = "grey80", alpha = 0.4, linetype = 3) + 
  geom_text(data = cluster_nested, 
            aes(x = cent_long, y = cent_lat, label = cluster), inherit.aes = FALSE) +
  geom_glyph(height = 2, width = 4) + 
  theme_void()

p2 <- ggplot() + 
  geom_sf(data = state_map, inherit.aes = FALSE,
          color = "grey80", alpha = 0.4, linetype = 3) + 
  geom_point(data = station_nested, aes(x = long, y = lat), size = 0.5) +
  ggforce::geom_mark_hull(data = cluster_nested |> tidyr::unnest(hull),
                          expand = 0, radius = 0,
                          aes(x = long, y = lat, group = cluster)) + 
  theme_void()

(p1 | p2) + plot_annotation(tag_levels = 'A')

Profile of aggregated precipitation at 639 weather stations in Australia. Subplot A shows the glyph map of the weekly averaged precipitation of each cluster. The group number is printed in the middle of the y minor axis and can be used as a reference line to read the magnitude. Subplot B shows the station membership of each cluster.

## clean <-  climate_full |>
##   filter(!id %in% c("ASN00067033", "ASN00072091", "ASN00059040",
##                     "ASN00097053", "ASN00040856", "ASN00015548")) |>
##   face_temporal() |>
##   mutate(month = lubridate::month(date)) |>
##   group_by(month) |>
##   summarise(tmax = mean(tmax, na.rm = TRUE),
##             tmin = mean(tmin, na.rm = TRUE),
##             diff = mean(tmax - tmin, na.rm = TRUE)) |>
##   ungroup(month) |>
##   face_spatial() |>
##   mutate(temp_diff_var = ts |> pull(diff) |> var(na.rm = TRUE)) |>
##   face_temporal() |>
##   mutate(dummy_date = as.Date(glue::glue("2021-{month}-01"))) |>
##   unfold(name)
## 
## df_id <- clean$id |> unique()
## p <- purrr::map(1:length(df_id), function(i){
##   dt <- clean |> filter(id == df_id[i])
##   dt |>
##     ggplot(aes(x = dummy_date, ymin = tmin, ymax = tmax, group = id)) +
##     geom_ribbon(color = "grey50", fill = "grey50", alpha = 0.5) +
##     theme_bw() +
##     scale_x_date(date_labels = "%b", date_breaks = "1 month") +
##     labs(x = "Month", y = "Temperature") +
##     ylim(-10, 50) +
##     ggtitle(paste0(stringr::str_to_title(unique(dt$name))))
## })
## 
## nested <- face_spatial(clean)
## domain <- nested$temp_diff_var
## pal <- colorNumeric(
##   colorspace::sequential_hcl(
##     "Rocket",  n = 7, cmax = 90, rev = TRUE,
##     c2 = 40, l2= 85, c1 = 20, l1 = 30),
##   domain = domain)
## 
## leaflet(nested) |>
##   addTiles() |>
##   addCircleMarkers(color = ~pal(domain), radius = 0.1,
##                    fillOpacity = 1, opacity = 1,
##                    popup = popupGraph(p, width = 300, height =200))

Same as Figure 11 except the temperature variation is now shown as a popup in the leaflet map.

## library(tidyverse)
## library(cubble)
## all_stations <- rnoaa::ghcnd_stations() |>
##   filter(str_starts(id, "ASN")) |> # Australian stations start wiht "ASN"
##   filter(last_year >= 2020) |>
##   mutate(wmo_id = as.numeric(wmo_id), name = str_to_lower(name)) |>
##   select(-state, -gsn_flag) |>
##   select(id, longitude, latitude, elevation, name,
##          wmo_id, element, first_year, last_year) |>
##   rename(long = longitude, lat = latitude, elev = elevation)
## 
## tmax_stations <- all_stations |>
##   filter(element == "TMAX", first_year < 1970, !is.na(wmo_id))
## 
## raw_tmax <- all_stations |>
##   rowwise() |>
##   mutate(ts = list(rnoaa::meteo_pull_monitors(
##     monitors = id, var = "TMAX",
##     date_min = glue::glue("{first_year}-01-01"),
##     date_max = glue::glue("{last_year}-12-31")
##     ) |>
##       select(-id)
##     )
##   )
## 
## historical_tmax <- raw_tmax |>
##   select(-element) |>
##   unnest(ts) |>
##   mutate(tmax = tmax/10) |>
##   filter(lubridate::year(date) %in% c(1971: 1975, 2016:2020)) |>
##   as_cubble(index = date, key = id, coords = c(long, lat))
## 
## save(historical_tmax, file = here::here("data/historical_tmax.rda"))

## aus_stations <- all_stations |>
##   filter(element %in% c("PRCP", "TMAX", "TMIN")) |>
##   nest(element: last_year) |>
##   rowwise() |>
##   filter(nrow(data) == 3) |>
##   select(-data)
## 
## aus_climate_raw <- aus_stations |>
##   rowwise() |>
##   mutate(ts = list(
##     rnoaa::meteo_pull_monitors(
##       monitors = id, var = c("PRCP", "TMAX", "TMIN"),
##       date_min = "2016-01-01", date_max = "2020-12-31"
##       ) |>
##       select(-id)
##     )
##   )
## 
## climate_full <- aus_climate_raw |>
##   unnest(ts) |>
##   mutate(tmax = tmax/10, tmin = tmin/10) |>
##   cubble::as_cubble(key = id, index = date, coords = c(long, lat))
## 
## save(climate_full, file = here::here("data/climate_full.rda"))